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Quantized vortex states of weakly interacting Bose-Einstein condensate of atoms with attractive 
interatomic interaction in an axially symmetric harmonic oscillator trap are investigated using the 
numerical solution of the time-dependent Gross-Pitaevskii (GP) equation obtained by the semi- 
implicit Crank-Nicholson method. Collapse of the condensate is studied in the presence of deformed 
traps with a larger frequency along the radial as well as along the axial directions. The critical 
number of atoms for collapse is calculated as a function of vortex quantum L. The critical number 
increases with angular momentum L of the vortex state but tends to saturate for large L. 
PACS Number(s): 02.70.-c,02.60.Lj,03.75.Fi 



I. INTRODUCTION 

Recent experiments |y,|| of Bose-Einstein condensates 
(BEC) in dilute bosonic atoms employing magnetic traps 
at ultra-low temperatures have intensified theoretical in- 
vestigations on various aspects of the condensate || [f|. 
The properties of the condensate are usually described 
by the nonlinear mean-field Gross-Pitaevskii (GP) equa- 
tion J8j, which properly incorporates the trap potential 
as well as the interaction among the atoms. 

Two interesting features of BEC are (a) the collapse 
in the case of attractive atomic interaction g^ and (b) 
the possibility of the fornration of a vortex state in har- 
monic traps with cylindrical |9|-|l3| as well as spherical 
J14J symmetry. 

For attractive interatomic interaction [g[7), the con- 
densate is stable for a critical maximum number of atoms. 
When the number of atoms increases beyond this crit- 
ical value, due to interatomic attraction the radius of 
BEC tends to zero and the maximum density of the con- 
densate tends to infinity. Consequently, the condensate 
collapses emitting atonrs until the number of atoms is re- 
duced below the critical number and a stable configura- 
tion is reached. With a supply of atonrs from an external 
source the condensate can grow again and thus a series 
of collapses can take place, which was observed experi- 
mentally in the BEC of 7 Li with attractive interaction 
Q. Theoretical analyses based on the GP equation also 
confirm the collapse 0. 

The study of superfluid properties of BEC is of great 
interest to both theoreticians PHTTJ] and experimentalists 
fl8| , [l9| . Quantized vortex state in BEC is intimately con- 
nected to the existence of superfluidity. Such quantized 
vortices are expected in superfluid Hejp However, due 
to strong interaction between the helium atoms there is 
no reliable mean-field description. On the other hand, a 
weakly interacting trapped BEC is well-described by the 
mean-field GP equation which is known to admit vor- 
tex solutions for a trap with cylindrical symmetry P,E6[ , 
that can be studied numerically. This allows for a con- 
trolled theoretical study of quantized vortices in BEC in 
contrast to superfluid Hejj. 



Many different techniques for creating vortex states in 
BEC have been suggested JlJ], e.g., stirring the BEC by 
an external laser at a rate exceeding a critical angular 
velocity to create a singly quantized vortex line along 
the axis of rotation [pT|, spontaneous vortex formation 
in evaporative cooling p0[ , controlled excitation to an 
excited state of atoms |21| , and rotation of an axially 
symmetric trap [J22|. Moreover, quantized vortex states 
in BEC have been observed experimentally in coupled 
BEC comprising of two spin states of 87 Rb in spherical 
trap, where angular momentum is generated by a con- 
trolled excitation of the atoms between the two states 
|r9fl . Vortices have also been detected in a single-state 
BEC of 87 Rb in cylindrical trap, where angular momen- 
tum is generated by a stirring laser beam |18J . After the 
possibility of continuously changing the interaction be- 
tween cold 85 Rb atoms by a magnetic-field-induced Fes- 
hbach resonance |23j24|, one could experimentally form 
vortex states in repulsive condensates and study their col- 
lapse after transforming them to attractive condensates 
by such a resonance. Because of the intrinsic interest 
in BEC of vortex states in axially symmetric traps, in 
this work the formation of such a BEC is studied using 
the numerical solution of the time-dependent GP equa- 
tion with special attention to its collapse for attractive 
interatomic interaction. 

In general, a vortex line in a nonrotating trapped BEC 
is expected to be nonstationary. However, it is possible 
to have dynamically stable vortex BEC states in a non- 
rotating trap with low quanta of rotational excitation 
or angular momentum L per particle |||l5|,[l6],^2| . Vor- 
tex BEC states for large repulsive condensates with high 
quanta of rotational excitation are expected to be unsta- 
ble and decay to vortices with low quanta |11 l^ , p^| . [l7| . 
In the absence of vortex, the stable condensate in an axi- 
ally symmetric trap has a cylindrical shape. Such a BEC 
has the largest density on the axis of the trap. For purely 
attractive interaction, with the increase of the number of 
atoms the central density of this condensate increases 
rapidly leading to instability and collapse M. 

In the presence of vortex motion the region of largest 
density of the BEC with nonzero L is pushed away from 



the central axial region and the atoms have more space to 
stabilize. The vortex state of the condensate in a cylin- 
drical trap has the shape of a hollow cylinder with zero 
density on the axis of symmetry. Because of larger espa- 
cial extension of such a condensate, it can accommodate 
a larger critical number of atoms before the density in- 
creases too high to lead to collapse |J. The higher the 
angular momentum L in a BEC, the larger is the critical 
number of atoms. However, the increase of this critical 
number with L slows down as L increases. 

The present study is performed with the direct nu- 
merical solution of the time-dependent GP equation with 
an axially symmetric trap. In the time-evolution of the 
GP equation the radial and axial variables are dealt with 
in two independent steps. In each step the GP equa- 
tion is solved by discretization with the Crank-Nicholson 
rule complimented by the known boundary conditions 
p5[ . We find that this time-dependent approach leads to 
good convergence. There are several other iterative ap- 
proaches to the numerical solution of the time-dependent 
and time-independent GP equation for axially symmetric 
P,p| p0|j2^ , |27[ | as well as spherically symmetric || traps. 
Of the time-dependent methods, the approach of Refs. 
pi| uses alternative iterations in radial and axial direc- 
tions as in this study, whereas Ref. [g6| does not give the 
details of the numerical method employed and Ref. |27| 1 
employs a completely different scheme, e.g., uses alterna- 
tive iterations for the real and imaginary parts of the GP 
equation. However, Refs. P,fO[ do not provide enough 
details of the numerical scheme. Because of these a mean- 
ingful comparison of the present method with those of 
Refs. [§l|@,||l is not possible. 

In Sec. II we describe the time-dependent form of the 
GP equation including the vortex states for attractive in- 
teraction. In Sec. Ill we describe the numerical method 
for solving the time-dependent GP equation in some de- 
tail. In Sec. IV we report the numerical results for the 
collapse of the BEC with vortex quantum for attractive 
interaction and finally, in Sec. V we give a summary of 
our investigation. 



II. NONLINEAR GROSS-PITAEVSKII 
EQUATION 

At zero temperature, the time-dependent Bose- 
Einstein condensate wave function ^(r, r) at position r 
and time r may be described by the self-consistent mean- 
field nonlinear GP equation g. In the presence of a 
magnetic trap of cylindrical symmetry this equation is 
written as 



2m 



V 2 + V(r)+gN\y(r,T) 



ifi 



d_ 



*(r,r) =0. 



(2.1) 



Here m is the mass of a single bosonic atom, N the 
number of atoms in the condensate, V(r) the attractive 
harmonic-oscillator trap potential with cylindrical sym- 
metry, g — A^nfi a/m the strength of interatomic interac- 
tion, with a the atomic scattering length. A positive a 
corresponds to a repulsive interaction and a negative a to 
an attractive interaction. The normalization condition of 
the wave function is 



Jdr\V(r,t)\ 2 = l. 



(2.2) 



The trap potential with cylindrical symmetry may be 
written as V(r) = ^mu 2 (r 2 + X 2 z 2 ) where u is the an- 
gular frequency of the potential in the radial direction r 
and A is the ratio of the axial to radial frequencies. We 
are using the cylindrical coordinate system r = (r, 9, z) 
and in case of cylindrical symmetry the wave function 
is taken to be independent of 9 in the absence of vortex 
states of the condensate: 



\l/(r, r) = ip(r, z, t). 



(2.3) 



The GP equation with a cylindrically symmetric trap 
can easily accommodate quantized vortex states with ro- 
tational motion of the condensate around the z axis with- 
out any added complication. In such a vortex the atoms 
flow with tangential velocity LTi/(mr) such that each 
atom has quantized angular momentum Lh along z axis. 
This corresponds to an angular dependence of 



^(r, r) = ip(r, z, r) exp(iL9) 



(2.4) 



of the wave function, where exp(iL9) are the circular 
harmonics in two dimensions. Equation (2.3) is the zero 
angular momentum vers ion of (2.4). 



Substituting Eq. (2^4) into Eq. (2T), one obtains the 
following GP equation in partial-wave form with quan- 
tized angular momentum L along the z axis 
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dr 



ip{r,z,r) = 0, 



(2.5) 



with L = 0, 1, 2, ... The nonzero values of L c orr esponds 
to vortex states. The L 2 /r 2 term in Eq. (2.5) is the 
vortex contribution to the Hamiltonian of the GP equa- 
tion. This is also the centrifugal barrier term in the 
partial-wave linear Schrodinger equation. The limitation 
to cylindrical symmetry reduces the GP equation in three 
space dimensions to a two-dimensional partial differen- 
tial equation. We shall study numerically this equation 
in this paper to understand the effect of the L 2 /r 2 term 
to collapse in the case of attractive atomic interaction. 

It is convenient to use dimensionless variables defined 
by x = \/2r/l, y = \/2z/l, t = tuj, and 



0{x,y,t) 



<p(x,y,t) 
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ip(r,z,T), (2.6) 



where I = y/h/(muj). Although <fi(x,y,t) is the dimen- 
sionless wave function, for calculational purpose we shall 
be using ip(x, y,t) in the following. In terms of these 
variables Eq. (2.5) becomes 



d 2 Id d 2 L 2 1 / , - , 4 



dx 2 



x dx 
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dy 2 
<p(x,y,t) 



X 
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<f(x,y,t) = 0, (2.7) 



where n = Na/l. A reduced number of p art icles is defined 
|n|. The normalization condition (J2.2J) of the wave 



as 



function become 
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dy\<p(x,y,t)\ 2 x- 



(2.8) 



However, physically it could be more interesting to de- 
fine the reduced number of particles in terms of a geo- 
metrically averaged frequency luo — A 1 ' 3 ^ and a length 
Iq = ^Jfi/mujQ, so that a new reduced number fc(A) is 
defined via J2S 



k(X) 



N\a\ 

In 



iX 1 ^. 



(2.9) 



We shall study this number in the present paper. 

For a stationary solution the time dependence of the 
wave function is given by <p(x,y,t) = exp(—ifit)tp(x,y) 
where /i is the chemical potential of the condensate in 
units of T im. If we use this form of the wave function in 
Eq. (2.7), we obtain the following stationary nonlinear 
time-independent GP equation pi: 



d 2 Id d 2 L 2 1 / „ , 9 , 4 
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dx 2 



x dx 
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dy 2 
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<p{x,y) = 0. (2.10) 



Equation (2.1C) is the stationary vers ion of the time- 
dependent Eq. (2/7). However, Eq. (|2.7J) is equally useful 
for obtaining a stationary solution with trivial time de- 
pendence as well as for studying evolution processes with 
expli cit t ime dependence and we shall be directly solving 
Eq. (2/7) numerically in this paper. 

Two interesting properties of the condensate wave 
function are the mean-square sizes in the radial and axial 
directions defined, respectively, by 



= 2tt 



dx 



dyx\<p(x,y,t)\* 



(2.11) 



III. NUMERICAL METHOD 

To solve the time-independent GP equation we need 
the boundary conditions of the wave function as x — > 
and co and \y\ — > oo. For a confined condensate, for a 
sufficiently large x and |y|, <p(x,y) must vanish asymp- 
totically. Hence the cubic nonlinear term can eventually 
be n eglect ed in the GP equation for large x and \y\ and 



2„,2 



2.10) becomes 
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dx 2 x dx dy 2 x 2 4 
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ip(x,y) = 0. 
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(3.1) 



This is the equation for the free oscillator with cylindrical 
symmetry in partial-wave form. The wave function for 
a general state of this oscillator and the corresponding 
energy are given, respectively, by |2£| 



and 




M=(l + |£|+n x )+ [n y + -)X, 



(3.3) 



with L = 0, ±1,±2, ..., n x = 0,2,4,..., and n y = 
0,1,2,... Here H ny is the usual Hermite polynomial, 
and -F|L|. rix is another polynomial defined recursively 
P8| , P9| , and A/" is the normalization. The first few of 
these polynomials are: Hq(£) = l,Hi(£) = 2^, ^(C) = 
(4C 2 - 2),F 3 (0 = C(8C 2 - 12),Fo,o(0 = 1,^1,0(0 = 

& ^2,0(0 - € 2 ,^ 2 (o = (i - e), Fs,o(o - c 3 , f 1)2 (o = 

^(^ 2 — 4) etc. Q. In this paper we shall be interested 
in angular momentum (vortex) excitation, opposed to 
radial excitation via n x or axial excitation via n y , of 
the following normalized ground state wave function for 







<p(x,y) 



with energy 
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x l+\L\ e -(x^+Xy')/4 
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(3.4) 



(3.5) 



and 



(y 2 ) = 2n dx dyx- l y 2 Mx,y,t)\ 2 . (2.12) 

JO J-oo 



Solution (3.4) of Eq. (3.1) is a good starting point for 
a iterative method for solving the time-dependent GP 
equation (2.7) for small values of nonlinearity n as in 
this paper. Alternatively, to solve the GP equation for 
large nonlinearity n, one may start with the Thomas- 
Fermi approximation for the wave function obtained by 



setting all the derivatives in the GP equation to zero || , 
which is a good approxim atio n for large nonlinearity. 

Next we consider Eq. (2.7) as x —* 0. The nonlinear 
term approaches a constant in this limit because of the 
regularity of the wave function at x — 0. Then one has 
the following condition 



¥>(o,y) = o, 



(3.6) 



as in the case of the harmonic oscillator wave function 
( |3.4|) . Both the small- and large- x behaviors of the wave 
function are necessary for a n ume rical solution of the 
time-dependent GP equation (2/7). The large- x and 
large- 1 y | behaviors of the wave function are given by Eq. 
(0, e-g-, 



observables using a finite difference scheme in this case 
does not lead to a tridiagonal set of equation but rather 
to a unmanageable set of equations |25|] . 

To circumvent this problem the full H operator in this 
case is conveniently broken up into radial and axial com- 
ponents H x and H y , respectively, where H x contains the 
terms dependent on x and H y the terms dependent on 
y with the nonlinear term 8\?2irn\(p(x,y)/x\ 2 involving 
both x and y contributing equally to both. Specifically, 
we take 



H x — — 



dx 2 
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x dx 



L 2 -l 



— + 4V2irn 
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V(x,y,t) 



x 
(3.12) 



lim (p(x, y) — > e x ' 4 , 
lim <p{x,y)^e- Xy2/4 . 



(3.7) 
(3.8) 



A convenient way to solve Eq. (2.7) numerically is to 
discretize it in both space and time and reduce it to a 
set of algebraic equations which could then be solved by 
using the known asymptotic boundary conditions. The 
method of solution using one space d eriv ative is well un- 
der control ||,^5[ . The GP equation (^7]) can be written 
formally as 



d 



(3.9) 



where H is the time-independent quantity in the square 
brackets of Eq. (2.7). The integration in time is ef- 
fected via the following semi-implicit Crank-Nicholson 
algorithm §§ 



<P" 



1 — io n 1 



(3.10) 



where A is the constant time step used to calculate the 
time derivative, ip n is the discretized wave function at 
time t n = nA, and where the space variables x and y are 
suppressed. The derivatives in the operator H are dis- 
cretized by the finite difference scheme |25| . The formal 
solution to Eq. ( 3.10 ) is given by 



^ = J-^/y 



1 + iAH/2 ' 



(3.11) 



so that if ip n is known at time t n one can find ip n+1 
at the next time step t n+ \. This procedure is used to 
solve the GP equation involving one space variable j|. 
In that case after proper disc retiza tion in space using a 
finite difference scheme Eq. (3.11) becomes a tridiago- 
nal set of equations in discrete space observables at time 
t n +i which is solved by Gaussian elimination method and 
back substitution |2q ] using the known boundary con- 
ditions (3.6), (|3.7[), and (|3.8[). Unfortunately, a similar 
straightforward discretization of Eq. (2.7) in two space 



H,. 



dy 2 



2, ,2 



y y 



aV2t, 



¥>{%,y,t) 



(3.13) 



with H — H x +H y . However, the numerical result of the 
present scheme is independent of a specific breakup. 

The procedure is then to define the unknown wave 
function on a two-dimensional mesh in the x — y plane. 
The time evolution is then performed in two steps. First 
the time evolution is effected using the operator H x set- 
ting H y — along lines of constant y with id(p/dt — H x ip. 
Next the time evolution is effected using the opera- 
tor H y setting H x = along lines of constant x with 
idtp/dt = Hyip. This procedure is repeated alternatively. 
This scheme is conveniently represented in terms of an 
auxiliary function ip n+ 2 by 



^4.1 = i-^i/y 



1 + iAHy/2 ■ 



n+1 l-iAH x /2 „ 



1 + iAH x /2 ■ 



so that 



<r 



n+1 _ 



(1 - iAH v /2) (1 - iAH x /2) 
(1+ iAH y /2) (l + iAH x /2) 



Y 



(3.14) 



(3.15) 



where n = 0, 1, 2, ... denotes the number of iterations. For 
a small time step A, if we neglect terms quadratic in A, 
Eq. (3.15) is equivalent to ( |3 .11 ). Hence for numerical 
purpose we have been able to reduce the GP equation 
in two space dimensions, x and y, into a series of GP 
equations in one space variable, either x or y. The GP 
equations in one space variable can be dealt with nu- 
merically in a standard fashion using Crank-Nicholson 
discretization and subsequent solution by the Gaussian 
elimination method. This scheme is stable independent 
of time step employed. 



The time-dependent GP equation (2.7) is solved by 
time iteration by mapping the solution on a two- 
dimensional grid of points N x x N y in x and y. First 



Eq. (2.7) with H x is discretized using the following fi- 
nite difference scheme along the x direction within the 
semi- implicit Crank-Nicholson rule [pTJ : 



i(tp 



n+l 
J,P 



^,p) 



1 
2h? 



Wi* - ^".p 1 + ^7-U 



<Pi 



+ (<Pj+l,p - 2t P),p T- ^j-l.j, 
, f/ n+l n+l \ i / n n \1 



4xjh 



2xi 



■2-y/27i 



j,pi 



«r+^j, (3.16) 



J,P 



where the discretized wave function tp^ p = (p(xj,y Pl t n ) 
refers to a hxed y = y p = ph,p = 1, 2, ..., iVj, at different 
x = Xj = jh,j = 1,2, ...,N X , and h is the space step. 
This proce dure results in a series of tridiagonal sets of 
equations (|3.16D in <p"+i p , ^J 1 , and ¥>"+i )P at time i n+1 
for each y p , which are solved by Gaussian elimination and 
back substitution pa] s tarting with the initial harmonic 
oscill ator solution (3.4) at to = and n = 0. Then 
Eq. (2.7) with H y is discretized using the following finite 
difference scheme along the y direction: 



i(tp 



n+l 
J,P 



<P?,p) 



A 



1 
2h? 



fcffl 



2w n+1 
rj,p 



+^;-i) 



+(4+i - 2 v1p + <P?*-i) 
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j.pi 



W 1 + <0> ( 3 - 17 ) 



where now </?!? refers to a fixed Xj = jh for all y p = ph. 
Using the solution obtained after x i terat ion as input, 
the discretized tridiagonal equations ( 3.17 ) along the y 
direction for constant x are solved similarly. This two- 
step procedure corresponds to a full iteration of the GP 
equation and the resultant solution corresponds to time 
t\ = A, and n = 1. This scheme is repeated for about 500 
times to yield the final s oluti on of the GP equation. The 
normalization condition ( |2.§| ) is preserved during time it- 
eration due to the unitarity of the time-evolution opera- 
tor. However, it is convenient to reenforce it numerically 
after each iteration in order to maintain a high level of 
precision. Also, the soluti on a t e ach time st ep w ill satisfy 
the boundary conditions (3.6), ( |3.7| ), and (3^8). At each 
iteration the strength of the non-linear term is increased 
by a small amount so that after about 500 time iterations 
the full strength is attained and the required solution of 
the GP equation obtained. The solution so obtained is 
iterated several times (between 20 to 50 times) until an 
equilibrated final result is obtained. This solution is the 
ground state of the condensate corresponding to the spe- 
cific nonlinear constant k and L. 

We found the convergence of the two-step iteration 
scheme to be fast for small \n\. However, the final con- 
vergence of the scheme breaks down if \n\ is too large. 
For an attractive interaction there is no such problem 
as the GP equation does not sustain a large nonlinearity 
\n\. Typical values of the parameters used in this paper 



for discretization along x and y directions are N x = 400, 
N y = 800, respectively, with x m ax = 8, |y|max = 8, 
and A = 0.05 for A > 0.5. For smaller A(< 0.5) the 
wave function extends further along the y axis and larger 
|ymax| and N y — 800 are to be employed for obtaining 
converged result. The above choice of parameters cor- 
responds to a typical space step of h — 0.02 along both 
radial and axial directions. These parameters were ob- 
tained after some experimentation and are found to lead 
to good convergence. 

As the time dependence of the stationary states is triv- 
ial — ip(x, y, t) — ip(x, y) exp(—ifj,t) — the chemical poten- 
tial /j, can be obtained from the propagation of the con- 
verged ground-state solution at two times, e.g., (p(x, y, t n ) 
and (p(x,y,t n + n /). From the numerically obtained ra- 
tio (p(x,y,t n )/(p(x,y,t n+n >) = exp(i^n'A), n can be ob- 
tained as the time step A is known. In the calculation 
of (i an average over relatively large values of n' leads to 
stable result. 



IV. NUMERICAL RESULT 

Using the numerical method described in Sec. Ill we 
present results for the numerical solution of the time- 
dependent GP equation in the following for attractive 
interatomic interaction with special attention to the col- 
lapse of the condensate. To assure that we are on the 
correct track, using the present program first we solved 
the GP equation for the spherically symmetric case with 
A = 1 and L = 0, and compare with the calculation of 
Rcf. |3(|. As an additional check we also solved the GP 
equation in two space di men sions with A = and without 
the d 2 /dy 2 term in Eq. (2.7) and compare with the calcu- 
lation of Ref . plf . In both cases the present calculation 
agrees with these previous ones. 

Before describing the results for nonzero L first we 
compare the present results for L = with those of Ref. 
pq| for a cylindrically symmetric trap. For the spher- 
ically symmetri c ca se A = 1, and the critical number 
k c (X) of Eq. (2.9) for collapse is found to be 0.575 
in agreement with Refs. |p|, p6|j30| . In a recent exper- 
iment using A = 0.3919, the critical reduced number 
for collapse for an attractive condensate of 85 Rb atoms 
formed using a Feshbach resonance was found to be 
k c = 0.459 ± 0.012 ± 0.054 Q. In their calculation 
Gammal et al |§ obtained k c = 0.550 for A = 0.3919. 
In the present calculation we obtain k c = 0.553 in excel- 
lent agreement with Ref. |26| using an entirely different 
numerical routine. However, the disagreement with the 
experimental result fl32] remains. We also calculated the 
critical number fc c (A) for some other values of A. For 
A = 5,2,0.2 we obtain k c = 0.50,0.56, and 0.52, respec- 
tively, compared to 0.498, 0.561, and 0.509 obtained in 
Ref. p6| . The small difference between the results of the 
two calculations seems to be a consequence of numerical 
error. Also, as in Ref. f[26l we note that for A not so 



different from unity (5 > A > 0.2) the critical reduced 
number for collapse k c (X) satisfies k c (X) rs fc c (l/A), and 
attains a maximum at A = 1 corresponding to the spher- 
ically symmetric situation. However, this symmetry is 
broken for large values of A, e.g., for A > 5 while we 
have k c (X) < fc c (l/A). Moreover, we find in the following 
that this symmetry is also broken for nonzero L, where, 
however, for A > 1 fc c (A) > k c (l/X). 

Next we comment on the discrepancy between the ex- 
perimental critical number of atoms for collapse for an 
attractive BEC of 85 Rb atoms formed using a Feshbach 
resonance |32| on one hand and the theoretical results of 
Ref. p6f and the present calculation on the other hand. 
In view of the success of the mean-field GP equation 
to explain many stationary results and time-evolution 
phenomena of the attractive BEC of 7 Li atoms with an 
almost spherical trap [|2||7|], it seems that this descrip- 
tion is perfectly appropriate for attractive condensates. 
Hence, we do not believe that a relatively small deviation 
from spherical symmetry as in the experiment of Ref. 
p2| would invalidate the applicability of the GP equa- 
tion to an attractive condensate. Whether the inclusion 
of higher order interaction terms in the mean-field GP 
equation could account for the observed data p(| yet re- 
mains to be established. To resolve the discrepancy we 
advocate further experimental study of collapse for at- 
tractive condensates after changing the trap symmetry 
(A). 

After the above preliminary comparative study, we 
present results for the numerical solution of the GP equa- 
tion ( |2.7| ) for nonzero L = 0, 1, 2, .., 8 and A = \/8 and 
1/a/8 for different fc(A). We recall that A = a/8 corre- 
sponds to the experiment of Ensher et al. M for the BEC 
of 87 Rb atoms. These two possibilities of A correspond 
to axial compression (A > 1) and elongation (A < 1) of 
the condensate, respectively. For each L we increase k 
from and calculate the chemical potential /j,. With the 
increase of k the wave function becomes more and more 
localized in space and beyond a certain value of k the 
density at the peak of the wave-function diverges and no 
stable normalizable solution of the GP equation with a 
well defined // can be obtained. 

In Fig. 1 we plot [i vs. k(X) for A = a/8 and 1/a/8 
for different L. We also exhibit the result for the spher- 
ically symmetric case A = 1 (L = 0) for comparison. 
The curves are plotted for all allowed values of k for 
the ground state in each case. The curves go up to a 
maximum critical value k c of k which defines the crit- 
ical number N c of atoms in that particular case via 
k c = N c \a\/lo. We find that (i) k c for a particular A in- 
creases with L and that (ii) k c for a particular nonzero L 
increases as A increases from l/v8 to v8, which demon- 
strates the breakdown of the numerically noted symme- 
try fc c (A) « fc c (l/A) for L = 0. To demonstrate these 
two effects in an explicit fashion we plot in Fig. 2 k c 
vs. L for A = \/8, 1, and 1/a/8. The three curves in- 
tersect approximately at L — which demonstrates that 



k c {X = a/8) « k c (X = 1/a/8) < k c (X = 1) for L = 0, with 
fc c (A = a/8) = 0.54, k c (X = I/a/8") = 0.55, and k c (X = 
1) = 0.575. However, this symmetry is broken for non- 
zero L while k c {X = a/8) > k c {X = 1) > k c (X = 1/a/8). 
The critical number fc c (A) increases with L for all A, and 
we see from Fig. 2 that this rate of increase slows down 
as L increases. 

In Figs. 3 and 4 we plot the wave function \ 4>(x, y)\ = 
\tp(x,y)/x\ in dimensionless variables of Eq. ( |2.6| ). In 
Figs. 3 (a) — (c) we show the wave function for A = 1/a/8 
and L = 0, 2 and 4, respectively, where the parameter k is 
chosen to be very close to the critical value k c for collapse. 
The nature of the wave function is qualitatively different 
for zero and nonzero L. For L = the wave function is 
peaked on the y axis; whereas for nonzero L it is zero on 
the y axis and is peaked at some finite x. In all cases 
the peak is sharp and the density of atoms is very large 
on the peak. The BEC collapses with a slight increase in 
the parameter k. For smaller k the wave function has a 
much broader maximum. When k approaches k c a sharp 
maximum of the wave function appears very rapidly. To 
illustrate this in Fig. 3 (d) we plot the L = 2 wave 
function for k = 2.5. If we compare this with the wave 
function of Fig. 3 (b) for L — 2 and k = 2.58 sw k c , the 
change in the shape is explicit. 

In Figs. 4 (a) — (e) we plot the wave function for 
A = a/8 and for L = 0,2,4,6, and 8, respectively, for 
k rs k c . If we compare Figs. 3 with 4 for same L we 
find that for A = 1/a/8 the wave functions extend over 
a larger region along the y axis compared to those for 
A = a/8. This is apparent if we compare Fig. 3 (a) 
with 4 (a), and is expected as A = a/8 corresponds to a 
stronger harmonic oscillator potential in the y direction 
responsible for axial compression. From Figs. 3 and 4 
we find that for both A, the peak in the wave function 
moves further away from the y axis as L increases. 

To understand some aspects of the variation of k c with 
L and A exhibited in Fig. 2, we plot in Figs. 5 (a) and (b) 
the mean-square sizes (x 2 ) and (y 2 ) vs. k for different L 
and for A = 1/a/8 and a/8, respectively. The results for 
vortex states {L > 0) in the spherically symmetric case 
with A = 1, remain between the those for A = 1/a/8 and 
v8 and are not explicitly shown here. For nonzero L the 
system acquires a positive rotational energy L 2 jx 2 which 
allows it to move away from the axial direction y. For L = 
the region of highest density is the y axis. For L^0 the 
density is zero on the y axis and has a maximum at some 
finite x. Consequently, the condensate has the shape of a 
hollow cylinder. Because of vortex motion the condensate 
swells and has more space to stabilize. Hence for L > 
the density does not go to an unbearable level with the 
same number of atoms as for L = 0, and k c increases with 
L for all A. However, for all L and A with the inc rease 
of nonlincarity k (or n) in the GP equation (2.7), the 
attractive nonlinear interaction term takes control and 
eventually the mean-square sizes ({x 2 } and (y 2 )) reduce 
as can be seen from Figs. 5. This eventual shrinking 



in size with the increase of the number of atoms for all 
L and A together with the outward push due to vortex 
motion for nonzero L takes the density of the BEC at 
the maximum of the wave function to an unbearably high 
level at some critical value k c of k leading to collapse. 

Although, for a fixed A, k c increases with L, the rate 
of increase slows down for large L. As k (or n) increases 
sufficiently for large L (> 8), the nonlinear term contain- 
ing n becomes the deciding factor in the GP equation 
and the L 2 /r 2 term starts to play a secondary role. Con- 
sequently, the increase in the critical number k c with L 
slows down as L increases and the number k c tends to 
saturate as can be seen clearly in Fig. 2. In all cases 
(A = V8) 1 and l/v8) this tendency of saturation is 
visible beyond L — 4. 



V. SUMMARY 

In this paper we present a numerical study of the time- 
dependent Gross-Pitaevskii equation under the action 
of a harmonic oscillator trap with cylindrical symme- 
try with attractive interparticle interaction to obtain in- 
sight into the collapse of vortex states of BEC. The time- 
dependent GP equation is solved iteratively by discretiza- 
tion using a two-step Crank-Ni cho lso n sc heme. We ob- 
tain the boundary conditions (3.6), (3.7), and (|3.8[) of 
the solution of the dimensionless GP equation ( |2.7D and 
use them for its solution. The solution procedure is ap- 
plicable for both attractive and repulsive atomic interac- 
tions as well as for both stationary and time evolution 
problems. It is expected that numerical difficulty should 
appear for large nonlinearity or large values of reduced 
number of particles k and large vortex quantum L. For 
medium nonlinearity, as in this paper, the accuracy of the 
time-independent method can be increased by reducing 
the space step used in discretization. 

The ground-state wave function for each L is found to 
be sharply peaked for attractive interatomic interaction 
with the parameters set close to those for collapse. In 
the case of an attractive interaction, the mean square 
sizes (x 2 ) and (y 2 ) decrease as the number of particles 
in the condensate increases towards the critical number 
for collapse. Consequently, the density increases rapidly 
signaling the on-set of collapse beyond a critical reduced 
number k c . 

The presence of the quantized vortex states increases 
the stability of the BEC with attractive interaction. The 
critical number fc c (A) for L = is largest in the spheri- 
cally symmetric case A = 1. For vortex states (L ^ 0), 
k c (X) increases with A. As the vortex quantum L in- 
creases, k c also increases. However, in the present calcu- 
lation a tendency of saturation in the value of k c is noted 
with the increase of L. As the parameter n or k in the 
GP equation increases, the nonlinear term starts to play 
the dominating role in the GP equation compared to the 
angular momentum term L 2 /x 2 . Once this happens, the 



rate of increase of k c with L slows down and it is not 
unlikely that the critical number attains a limiting max- 
imum value for a larger L(> 8) than those considered in 
this paper. This and other investigations on the collapse 
of vortex states are welcome in the future. 
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Figure Caption: 

1. Chemical potential /i vs. reduced number k for dif- 
ferent A and L. The curves are labeled by their respective 
L value. 

2. Critical reduced number k c vs. L for A = y/8 (full 
line with x), 1 (dashed-dotted line with +) and 1/V& 
(dashed line with *) . The lines are polynomial fits to the 
points. 

3. The wave function \<p(x, y)\ = \ip(x,y)/x\ vs. x and 
y for A = l/v/8 and for (a) L = 0, k = 0.54 (b) L = 2, 
k = 2.58 (c) L = 4, k = 4.00 and (d) L = 2, k = 2.50. 

4. Same as figure 3 for A = \/8 and for (a) L = 0, 
k = 0.53, (b) L = 2, k = 3.15, (c) L = 4, k = 5.13, (d) 
L = 6, k = 6.71, and (e) L = 8, k = 8.12. 

5. Mean square sizes (x 2 ) (full line) and (y 2 ) (dashed 
line) vs. reduced number k for (a) A = l/v8 and (b) 
A = V8. 
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3. The wave function \<p(x,y)\ vs. x and y for A = 
I/1/8 and for (a) L = 0, k = 0.54 (b) L = 2, k = 2.58 (c) 
L = 4, k = 4.00 and (d) L = 2, k = 2.50. 






4. The wave function \(f>{x,y)\ vs. x and y for A = \/8 
and for (a) L = 0, fc = 0.53, (b) L = 2, fc = 3.15, (c) 
L = 4, fc = 5.13, (d) L = 6, k = 6.71, and (d) L = 8, 
fc = 8.12. 
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